Homework 2: Secant method

This assignment teaches implementing simple algorithms with loops and conditionals.

Define a function secant that takes as input

a function f
a lower bound xlo
an upper bound xhi
a resolution dx

that finds a root of f in the interval [xlo:xhi] with a resolution of dx. That is, the function secant should iteratively reduce the size of the interval until it is smaller than dx, bracketing a root of f.

Note:

Check initially whether there is actually a root inside the interval, i.e. whether f(xlo) and f(xhi) have different sign
Use the secant method; alternatively, use the bisection method which is slower, but slightly easier to implement
If you pass a function as argument to a function, it is customary to not specify a type constraint for this argument, as in function secant(f, xlo::Real, xhi::Real, dx::Real)
To test performance (in the sense of number of iterations required), use the info command to output a message during each iteration
As a test, find the value of sqrt(2) to 100 digits. Hint: Solve x^2-2==0, and use the BigFloat type.



In [1]:

    
typeof(+)









    Out[1]:





Function



In [10]:

    
midpoint(xlo::Real,xhi::Real)=0.5*(xlo+xhi)









    Out[10]:





midpoint (generic function with 1 method)



In [29]:

    
function secant(f,
    xlo::Real,xhi::Real,
    dx::Real)
    # check various conditions
    if xlo == xhi
        throw(DomainError())
    end
    if xlo > xhi
        return secant(f,xhi,xlo,dx)
    end
    flo,fhi = f(xlo),f(xhi)
    if sign(flo*fhi) > 0
        throw(DomainError())
    end
    if abs(xlo) <= dx return xlo end
    if abs(xhi) <= dx return xhi end
    # Initial two guesses based on midpoint method. 
    x0=midpoint(xlo,xhi)
    f0=f(x0)
    info("x=$x0 | f=$f0")
    if abs(f0) <= dx return x0 end
    if f0 < 0
        x1 = x0 < 0 ? midpoint(xlo,x0) : midpoint(x0,xhi)
    else
        x1 = x0 < 0 ? midpoint(x0,xhi) : midpoint(xlo,x0)
    end
    info("x=$x1 | f=$(f(x1))")
    if abs(f(x1)) <= dx return x1 end
    # and now the loop!
    while abs(f(x1)) > dx
        x0,x1 = x1,x1 - f(x1)*(x1-x0)/(f(x1)-f(x0))
        info("x=$x1 | f=$(f(x1))")
    end
    return x1
end









    Out[29]:





secant (generic function with 1 method)



In [35]:

    
with_bigfloat_precision(400) do
    secant((x)->x^2-2,BigFloat(0.5),BigFloat(2.),BigFloat(10.0^(-100)))
end









    



INFO: x=1.2500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 | f=-4.3750000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-01
INFO: x=1.6250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 | f=6.4062500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-01






    Out[35]:





1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350175239968496293396






    



INFO: x=1.4021739130434782608695652173913043478260869565217391304347826086956521739130434782608695652173913043478260869565217391303 | f=-3.3908317580340264650283553875236294896030245746691871455576559546313799621928166351606805293005671077504725897920604915751e-02
INFO: x=1.4133752244165170556552962298025134649910233393177737881508078994614003590664272890484739676840215439856373429084380610416 | f=-2.3704750055600501532639911812769743012870307398251082195268961382631370286447337461200519582657800669784592375801372437267e-03
INFO: x=1.4142171472137510396451344607707235930135846964236207374549487108400332686443027446631549764347102855558635985583587468807 | f=1.0139473400379956141406034488125659614053289014083987416575707725382882635334107224124240447043756906468024278024531230222e-05
INFO: x=1.4142135613102445402369252307717469547019258230346008829978111759964780573161412860528082303267007539911490085405847082362 | f=-3.0061952068452086966522608034736850536047173398361598168958260150619713032195808524714453767134084472195930555803966260185e-09
INFO: x=1.4142135623730937017120453891539138048225234649996533611509123887641140944275433892741760891922413543633118168461394280963 | f=-3.8101448866735407380831360860691007133762125252173046218399459530735453512230012503214986350769986019329937848957840445516e-15
INFO: x=1.4142135623730950488016892304115020820885762754880409766226633997792402545888716828391961774192887418224352651599184572156 | f=1.4317549130390074341884289325222480437920034419633290444562554569430320975725264864230976742275066687868397548656188023595e-24
INFO: x=1.4142135623730950488016887242096980785694307876000544567699984040296199470577871771324796742709608064162388566794955193654 | f=-6.8189920761066241888012871296255599892876719662803027884739871414234053044880737776171499088669598095755371321943300854635e-40
INFO: x=1.4142135623730950488016887242096980785696718753769480731766797379475851523319711158261102618135177262502102474282952659514 | f=-1.2203906758674650598501283064687850372236992054062166595072768940005831675790196324201868176667745690680478368442215109566e-64
INFO: x=1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350175239968496293396 | f=1.04022929356183148561185375120515176961354446918268928571109431669636469496537946540804333365334740504478630712730151375e-104



In [ ]: